
c = 299792458.d
Sigma = 3.d-4

;; Load in OH5 dust, Spitzer response function, and the {Draine:2007} j_\nu.
bandpass = rd_tfile('tpeb_analysis/IRAC_ch4trans_full.txt',2,-1,/convert)
sband = bandpass[1,*]
irac_lam = bandpass[0,*]

;; Interpolate to a common (frequency) grid.
nu_range = c / [11.d-6, 6.d-6]

dnu = 1.d11

nu = dindgen(230)*dnu + 2.7d13
print,m4_stat(nu)

print,'Delta nu = ',dnu,'  Delta lambda = ',-c/3.8075d13/3.8075d13*1.d11


;; Interpolate now...
s_band = interpol(sband, c/(irac_lam*1.d-6), nu)

plot,nu,s_band,psym=2


plots,c/6.d-6,0.35,psym=2,thick=5,color=cgColor('Red')
plots,c/7.d-6,0.35,psym=2,thick=5,color=cgColor('Red')
plots,c/8.d-6,0.35,psym=2,thick=5,color=cgColor('Red')
plots,c/9.d-6,0.35,psym=2,thick=5,color=cgColor('Red')
plots,c/10.d-6,0.35,psym=2,thick=5,color=cgColor('Red')
plots,c/11.d-6,0.35,psym=2,thick=5,color=cgColor('Red')



openw,lun,'/d1/BGPS/bgps_svn/papers/IRDC_prior/tables/kband.tex',/GET_LUN
;;==========================================================================
;; Table Header
printf,lun,'%% KBAND Table'
printf,lun,''
printf,lun,'\begin{deluxetable}{lcccccc}'
printf,lun,'  \tablecolumns{7}'
printf,lun,'  \tablewidth{0pc}'
printf,lun,'  \tabletypesize{\small}'
printf,lun,'  \tablecaption{IRAC Band 4 $\langle\kappa\rangle_{_\mathrm{band}}$\label{table:kband}}'
printf,lun,'  \tablehead{'
printf,lun,'    \multicolumn{2}{c}{Dust Emission Model\tablenotemark{a}} & \colhead{} &\multicolumn{4}{c}{$\kband$} \\'
printf,lun,'    \cline{1-2} \cline{4-7}'
printf,lun,'    \colhead{$q_{_\mathrm{PAH}}$} & \colhead{$U_{min}$} & \colhead{} & \colhead{OH5\tablenotemark{b}} & \colhead{D03-3.1\tablenotemark{c}} & \colhead{D03-5.5\tablenotemark{d}}  & \colhead{WD01-5.5\tablenotemark{e}}\\'
printf,lun,'    \colhead{(\%)} & \colhead{} & \colhead{} & \multicolumn{4}{c}{(cm$^2$ g$^{-1}$)}'
printf,lun,'  }'
printf,lun,'  \startdata'










;; Loop over DL07 dust models

umin = ['0.10','1.00','10.0']
modl = ['00','30','60']
qpah = ['0.47','2.50','4.58']

FOR j=0,n_elements(modl)-1 DO BEGIN
   FOR i=0,n_elements(umin)-1 DO BEGIN
      
      readcol,'./ancillary/U'+umin[i]+'_1e6_MW3.1_'+modl[j]+'.txt',$
              dl07_lam,jnu,format='f,x,f',skipline=61
      
      ;; Interpolate to a common (frequency) grid.
      j_nu = interpol(jnu, c/(dl07_lam*1.d-6), nu)
      
      kap_str = ''
      
      FOR kap=0,3 DO BEGIN
         
         CASE kap OF
            0: BEGIN
               OH5 = mrdfits('tpeb_analysis/OH5_opacities.fits',1,hdr)
               kappa = interpol(OH5.kappa2, c/(OH5.lambda2*1.d-6), nu)
               bgps_opac = 1.14d
            ENDCASE
            1: BEGIN
               readcol,'./ancillary/kext_albedo_WD_MW_3.1_60_D03.txt',$
                       lam,kabs,format='f,x,x,x,f',skipline=80
               kappa = interpol(kabs, c/(lam*1.d-6), nu)
               IF i+j EQ 0 THEN BEGIN
                  plot,lam,kabs,/xlog,/ylog,xr=[1,2000],/xst,thick=5,$
                       yr=[1d-2,1d4],ytickformat='exponent10'
                  oplot,OH5.lambda2,OH5.kappa2,thick=3,color=cgColor('Red')
                  vline,299792.458d/271.1d,color=cgColor('Cyan'),/xlog,/ylog
                  wait,0.05
               ENDIF
               bgps_opac = (interpol(kabs,lam,299792.458d/271.1d))[0]
           ENDCASE
            2: BEGIN
               readcol,'./ancillary/kext_albedo_WD_MW_5.5A_30_D03.txt',$
                       lam,kabs,format='f,x,x,x,f',skipline=80
               kappa = interpol(kabs, c/(lam*1.d-6), nu)
               IF i+j EQ 0 THEN BEGIN
                  plot,lam,kabs,/xlog,/ylog,xr=[1,2000],/xst,thick=5,$
                       yr=[1d-2,1d4],ytickformat='exponent10'
                  oplot,OH5.lambda2,OH5.kappa2,thick=3,color=cgColor('Red')
                  vline,299792.458d/271.1d,color=cgColor('Cyan'),/xlog,/ylog
                  wait,0.05
               ENDIF
               bgps_opac = (interpol(kabs,lam,299792.458d/271.1d))[0]
            ENDCASE
            3: BEGIN
               readcol,'./ancillary/kext_albedo_WD_MW_5.5B_30.txt',$
                       lam,kabs,format='f,x,x,x,f',skipline=43
               kappa = interpol(kabs, c/(lam*1.d-6), nu)
               IF i+j EQ 0 THEN BEGIN
                  plot,lam,kabs,/xlog,/ylog,xr=[1,2000],/xst,thick=5,$
                       yr=[1d-2,1d4],ytickformat='exponent10'
                  oplot,OH5.lambda2,OH5.kappa2,thick=3,color=cgColor('Red')
                  vline,299792.458d/271.1d,color=cgColor('Cyan'),/xlog,/ylog
                  wait,0.05
               ENDIF
               bgps_opac = (interpol(kabs,lam,299792.458d/271.1d))[0]
            ENDCASE
         ENDCASE
         
         ;; Calculate!
         
         I1 = total( 1./nu * s_band * j_nu * exp(-Sigma*kappa) * dnu )
         I0 = total( 1./nu * s_band * j_nu * dnu )
         
         ratio = I1 / I0
         ;; print,ratio
         
         
         ;; Spit out the answer.
         
         kband = -1./Sigma * alog(ratio)
         
         print,kband,bgps_opac,kband/bgps_opac
         
         kap_str += ' & '+string(round(kband),format="(I0)")
         
      ENDFOR
      
      ;; Print to table file
      sq = (i EQ 0) ? qpah[j]+'......' : ''
      
      printf,lun,sq+' & '+string(umin[i],format="(F0.1)")+' & '+kap_str+' \\'
      

   ENDFOR
   printf,lun,'\hline'
ENDFOR



printf,lun,'  \enddata'
printf,lun,'\tablenotetext{a}{From \citet{Draine:2007}}'
printf,lun,'\tablenotetext{b}{Extinction from \citet[][Table 1, Column 5]{Ossenkopf:1994}}'
printf,lun,'\tablenotetext{c}{$R_V = 3.1$ extinction from \citet{Draine:2003}}'
printf,lun,'\tablenotetext{d}{$R_V = 5.5$, extinction from \citet{Draine:2003}}'
printf,lun,'\tablenotetext{e}{$R_V = 5.5$, Model B, extinction from \citet{Weingartner:2001}}'
printf,lun,'\end{deluxetable}'


close,lun
free_lun,lun

END
